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Аннотация 

Введение. В области вычислительной математики известно множество способов аппроксимации модели 
механики жидкости. Учеными выработаны методы и оценки критериев качества аппроксимации, таких как 
устойчивость и сходимость. Комбинация подходов построения экономичных разностных схем, таких как 
расщепление по физическим процессам, регуляризация по Б. Н. Четверушкину, линейная комбинация 
разностной схемы «кабаре» и «крест» в совокупности ранее не реализовывалась и не оценивалась. Перед 
авторами стояла задача аппроксимировать каждую часть расщеплённой по физическим процессам модели 
гидродинамики наиболее адекватной схемой и далее исследовать корректность данного подхода. 

Материалы и методы. Математическая модель гидрофизических процессов замыкается эмпирическим 
уравнением состояния соленой воды. Выбираются значимые свойства, строится математическая модель. 
Разностные операторы аппроксимируют дифференциальные операторы. Строится алгоритм послойного 
моделирования переходных процессов. Алгоритм реализован в виде программы, которая, в основном, содержит 
поэлементные (массивно параллельные) операции. 

Результаты исследования. Получены математические модели гидродинамических процессов в водоемах, 
учитывающие три уравнения движения при наличии градиента плотности водной среды при отказе от 
гидростатического приближения. Апробирован новый способ вычисления поля давления с применением 
регуляризаторов по Б.Н. Четверушкину в уравнении неразрывности. Разработан программный модуль 
численного моделирования гидрофизических процессов движения воды с различной солёностью и плотностью. 
Это открытое программное обеспечение, допускающее не только переопределение эмпирических зависимостей 
(как алгебраических функций), но и подключение внешних моделирующих модулей для отображения 
зависимостей алгоритмически. 

Обсуждение и заключение. Разработанная модель гидрофизики, учитывающая свойства солёной воды и 
динамическую связь механического движения воды с солёностью, может применяться для изучения 
формирования неравновесного распределения параметров и идентификации наиболее стабильных параметров 
водной среды. Модель объясняет нисходящее движение кислорода, что позволит в будущем оценивать 
величины параметров водной среды, которые сложно измерить непосредственно. Она может быть использована 
в процедуре параметрической идентификации трудноизмеряемых параметров водной среды. 
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Введение. Одной из важных задач, связанных с экологией и безопасностью жизнедеятельности людей, 
живущих в прибрежных территориях, является прогнозирование и моделирование движения воды в морях и 
крупных региональных водоёмах. Кроме того, актуально изучение процесса переноса растворенных в водной 
среде веществ при учёте стратификации и зависимости плотности воды от многих переменных факторов. Такое 
прогнозное моделирование, вероятно, позволит не только оценивать качество вод, но и управлять им в 
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условиях климатических изменений и индустриальных воздействий. Более общая и главная цель — связать 
качество вод с численностью и видовым разнообразим естественных гидробиоценозов, обитающих в 
гидросфере. При решении этих задач необходимо учитывать гидродинамические характеристики водной среды, 
особенности внешних факторов воздействия, таких как неоднородность распределения температуры, соленость, 
насыщенность воды кислородом, количество растворенных в воде газов, кислотность. Эти параметры являются 
одними из ключевых параметров биологической активности водной экосистемы [1]. 

Акватория водоема также может рассматриваться как передаточная система, проводящая кислород от 
атмосферы к донным осадкам. Однако известны случаи образования химических градиентов большой 
величины на относительно небольших перепадах глубин — в пограничных слоях при относительно небольшой 
их величине за пределами этих зон. Причина появления таких участков при общей стратификации воды по 
плотности во всемирном тяготении с одной стороны и радиационном воздействии солнца на воду, приводящего 
к её нагреванию, с другой. Эти процессы могут привести к снижению скорости продукции, деструкции и 
рециклинга биогенных элементов и биоорганизмов вплоть до остановки этих процессов, а также 
предопределять биоразнообразие гидробионтов в целом и видовой состав в частности [2]. Температурная 
стратификация существенно воздействует на распределение организмов в толще воды, на перенос и осаждение 
вредных для биоорганизмов примесей. Повышение температуры поверхностных вод приводит к нарушению 
вертикального водообмена и, соответственно, к уменьшению аэрации глубоководной зоны, снижению 
растворимости и концентрации кислорода в воде. Стратификация по плотности, температуре и химическому 
составу ограничивает конвективный подъем в поверхностные слои воды биогенных элементов, углекислоты и 
продуктов неполного окисления органических веществ, поступающих в гиполимнион (холодные, солёные, 
плотные слои воды) в результате седиментации (оседания) сестона. С начала стратификации и до момента ее 
окончания поверхностный слой обедняется, а гиполимнион, напротив, обогащается этими веществами. В 
результате физико-химическая стратификация приводит к возникновению неравномерного распределения по 
глубине ряда биологически значимых веществ и является причиной самоорганизации сложной структуры 
экологических ниш [3]. Математическое моделирование механических, химических и биологических 
процессов, протекающих в водных экосистемах, чрезвычайно актуально и связано с проблемами экологии и 
безопасности жизнедеятельности населения прибрежных территорий. 

Среди выдающихся ученых, внесших значительный вклад в изучение гидрологии и океанологии, следует 
отметить В.П. Дымникова, который занимался исследованием климата и океанологии, моделированием 
атмосферы и океана. А.С. Монин и М.Ю. Белевич исследовали процессы кинематики водной среды, 
турбулентности и микроструктуры океана. Советский ученый А.С. Судольский занимался исследованием 
динамики вод и береговых процессов в различных водоемах применительно к решению задач проектирования, 
строительства и эксплуатации конкретных гидротехнических сооружений или рационального хозяйственного 
использования водоемов в целом [4]. По мнению В.И. Вернадского, одним из важнейших проявлений жизни 
является газовый обмен организмов с окружающей средой, главным образом — процессы дыхания, основанные 
на потреблении кислорода [5]. Изучением гидробиологических процессов водоемов занимались различные 
выдающиеся отечественные ученые, среди которых С.В. Бруевич, работы которого посвящены разработке 
аналитических методов исследований, формулировании основ гидро- и биогидрохимии. Г.Г. Матишов и 
В.Г. Ильичев активно изучают условия оптимальной эксплуатации водных ресурсов, занимаются разработкой 
моделей транспорта загрязняющих веществ в водоемах и исследованием оценки их воздействия на биоресурсы 
водной среды [6]. 

Для изучения влияния перечисленных процессов водной среды разрабатывается комплекс взаимосвязанных 
математических моделей, основанный на использовании точных прогностических моделей и программной 
реализации экономичных численных методов, который позволит детально исследовать кинематику процесса, 
причинно-следственные связи и состояние объекта моделирования. Существующие методы и средства 
прогнозного моделирования состояния водной среды с учетом ряда биотических и абиотических факторов, 
включая процессы распределения кислорода, углекислого газа, солей, основаны на общих научных подходах, 
упрощенных математических моделях, обладающих низкой адаптивностью, отсутствием возможности 
моделирования нелинейных динамических процессов, свойственных большинству водных экосистем, 
некорректным заданием границ расчетной области. В некоторых случаях исследования сопровождаются 
формальным определением граничных условий, дают достаточно грубые и приближенные результаты 
моделирования. 

При моделировании процесса переноса веществ, основанного на уравнениях адвекции-диффузии, 
необходима хорошая аппроксимация адвективных слагаемых, представляющих собой градиенты давления, 
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плотности массы и полной энергии, импульса движения. Применение стандартных разностных схем при 
завышенных оценках параметров подобия приводит к потере точности вычислений из-за возрастания 
погрешности аппроксимации и усилению ограничений на шаг по времени в связи с условием устойчивости 
разностной схемы. В работах А.И. Сухинова, А.Е. Чистякова и др. [7, 8] показано, как эффективно использовать 
линейную комбинацию разностных схемам «кабаре» и «крест» с оптимальными значениями весовых 
коэффициентов для аппроксимации уравнения переноса. Эффективность этих методов достигается 
оптимизацией погрешности аппроксимации дискретной моделью сплошной среды точного решения задачи 
переноса вещества с постоянной скоростью. Исследования показали, что такой подход распространяется и на 
модели гидродинамики с переменной (знакопеременной) скоростью без эффекта сеточной вязкости. Ещё одно 
положительное качество такого рода аппроксимаций состоит в том, что с их помощью можно моделировать 
сложные структуры течения, к примеру, вихревые. В настоящее время многие исследователи используют 
подобные схемы для моделирования турбулентных течений. Сотрудниками ведущих зарубежных научно- 
исследовательских организаций, таких как З‘апЮга ОшуегзИу, Ппрепа| СоШезе Гоп4доп и др., а также 
сотрудниками Института вычислительной математики РАН Е.М. Володиным, 
А.В. Глазуновым, А.С. Грицуном, Н.Г. Яковлевым и др. [9] изданы труды, в которых математическое 
моделирование климатических изменений, гидродинамических и атмосферных процессов и явлений 
осуществлено на основе вихреразрешающих схем. Ещё больше снизить требование к шагу по времени и 
повысить пространственное разрешение модели при ограниченной компьютерной памяти позволяет 
квазигидродинамическое приближение сплошной среды. На практике в систему уравнений Навье-Стокса и 
неразрывности добавляется малое слагаемое, пропорциональное второй производной по времени от функции 
плотности. Данный подход позволяет сгладить нефизические флуктуации плотности массы и импульса, а также 
полной энергии, переносимые по пространственной сетке быстрее скорости звука. 

Существующие универсальные пакеты прикладных программ (например, пакет программ «Магз3а», 
Экоинтегратор, программный комплекс СНАВ!ЗМА, комплекс ЗАГМО, комплекс программ СНТОМ, 
САКОПМАГ, пакеты моделирования различных процессов аэрогидродинамики, программные комплексы 
РНОЕМ!С$, ЕГОЕМТ, СА РУМАМ!С$ ТООГ) не учитывают некоторые свойства моделируемых сложных 
систем, снижая таким образом точность и оперативность моделирования. К таким свойствам относятся: 
пространственная неоднородность движения водной среды, вихревые структуры течений. В математических 
моделях и в алгоритмах их численной реализации не учитывается вероятность значительного изменения 
глубины, плотность водной среды, что может привести к неустойчивости полученных численных решений. По 
этой причине такие специализированные программные пакеты могут быть использованы для моделирования 
ограниченного разнообразия гидрофизических процессов водных систем. Большая часть известного 
специализированного программного обеспечения (АРАМ, САГЗОНС, СВепз1, ТАЗСЙо\, [$С-3, РАМАСНЕ, 
ВЕМ$АО, ОАМ-ГУ, ЭКОЛОГ, ПРИЗМА, УГПТЕСОМ), предназначенная для моделирования процесса 
распространения загрязняющих веществ, взаимодействия гидробионтов, преимущественно ориентирована на 
однопроцессорные вычислительные системы, представленные, в основном, персональными компьютерами. В 
таких системах масштабируются на параллельные системы только единичные составляющие модули этих 
систем (например, ЕСОЗ1М и МАОЗ1Р). На практике в систему уравнений Навье-Стокса и неразрывности 
добавляется малое слагаемое, пропорциональное второй производной по времени от функции плотности. 
Данный подход позволяет сгладить нефизические флуктуации плотности массы, импульса и полной энергии, 
переносимые по пространственной сетке быстрее скорости звука. 

Материалы и методы. Успешность разработки математической модели гидрофизических и 
гидробиологических процессов зависит от наличия и проработанности тестовых примеров и задач для изучения 
устойчиво наблюдаемых в морях явлений, таких как вертикальное перемешивание и перераспределение 
солёности и кислорода, галоклин и теромоклин. Для изучения этих явлений в работе используется модель 
гидродинамики, учитывающая баланс массовых сил и трансграничных потоков [1, 2]: 


Р(,рь) в, 
и 1Ур+Ь. 
ро |+ | - 4» (ть) + (,рь)+(У.п)-+ра , 
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где ру — плотность потока, 0/0 — скорость изменения плотности кинетической энергии, 
08/0: — скорость изменения плотности внутренней энергии, Т — тензор напряжений Т1и=-рп, Б — 
массовая сила, й = й(1,х) — плотность потока тепла, 4 — удельный приток тепла за счёт излучения. Из-за 


того, что эти явления обычно описываются вертикальным распределением параметров по глубине, то 
целесообразно получить упрощённую модель, допускающую оперативную идентификацию ненаблюдаемых 
параметров. Выделим цилиндрическую область У с основаниями на дне и поверхности воды площадью 
поперечного сечения 5. Спроецируем скорости, потоки и силы на вертикальное направление в предположении, 
что частные производные по х, у от параметров на горизонтальное направление равны 0 за пределами цилиндра, 
те. предположим горизонтальную однородность параметров водной среды. Запишем систему (1) В 
консервативной форме так, чтобы она позволила определить плотность массы ( Р ), механический импульс ( ру ) 


и плотность полной энергии (Р(К +=) ), К=у?/2. В цилиндре выделим бесконечно малый объём и 


предположим, что на каждый такой объём, составляющий У, действуют: сила реакции опоры (дна водобма), 
равная гидростатическому давлению, и сила, аналогичная силе трения, вызываемая вязкостью жидкости и 
переносом импульса, не равная нулю при вертикальных движениях жидкости. 

Если пренебрегать горизонтальными перемещениями жидкости и предположить, что существенны только 
вертикальные движения, и тот факт, что плотность водной среды существенно зависит от солёности, принять 
изложенные ранее упрощения и соглашения, то уравнения гидромеханики [1-3] в компактной форме 
представимы системой дифференциальных уравнений в частных производных [2]: 


2(ру) (ру?) __@р 
О Ох Ох 


д д 9(ру) № вот 
5: [Р(в+у" 2) +5 [+(р(в+т?/2)) |= о 2 р &(Т-Т.,),Т = Хе) ==/с, (1) 
& 0(у5) "| э 
+ = и— |, 
0 Ох Ох\о дх 


где р = Г(р,Т,5), == Г(р,у,Г) — эмпирическое уравнение состояния морской воды и уравнение, замыкающее 
систему по внутренней энергии соответственно; р — полное гидродинамическое давление; р — локальная 
плотность водной среды; у — проекция вектор-функции скорости на вертикаль (ось 7 направлена вверх от дна 
к поверхности); & —Й объемная плотность внутренней энергии; р — давление газа, заключённого в 
элементарном объеме между соседними слоями; РЕ — объемная плотность обобщённой силы (суммы сил), 
приложенных к элементарным объёмам жидкости, помимо давления; 5 = На) — концентрация соли; 5 — 
площадь поперечного сечения цилиндрической выделенной области, в которой предполагается наиболее 


интенсивное протекание процесса апвеллинга и транспорта соли; Т — абсолютная температура воды; Т. 


епу 
температура воды внешней по отношению к выделенному объёму; К — теплопроводность воды; 

$ =9,8м/ с? — ускорение силы тяжести; ци — коэффициент, характеризующий интенсивность переноса 
импульса из-за вязкости. 


Граничные условия, характеризующие свойство непротекания жидкости через породу, составляющую дно 
водного объекта, для системы (1) могут быть записаны в виде равенств: 


д ОТ д 
(уп) =0, 2 =0, —=0, ^^ =0. (2) 
Где п — вектор нормали, направленный внутрь расчётной области. 


Международный стандарт уравнения состояния морской воды [10] определяет плотность р морской воды 
как функцию солености 5, температуры Ти гидростатического давления р, которая имеет вид: 
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(3) 


где К(з, Т‚р) — средний модуль упругости; цифра 0 соответствует одной стандартной атмосфере (101 325 Па). 


При численном моделировании на основе (1)-(3), также востребованной как (3), установим 
алгебраическую связь давления с плотностью, температурой и солёностью. Квадратное уравнение имеет два 
корня: 


_ =Ур+(-А)-р(5.Т,р)+ А: рь(5,Т,р) 


Р!2 р, ь (4) 
р, = 2В-(р(5,Т, Р)-Ро(5,Т, р), 
р = (—4ВК,(5,Т, р) + А? -2А+1-р(5,Т, р)? + (8В-Ку(5,Т,р)-2А?-+2А). ру(5,Т,р): (5) 
р(5,Т, р) + (А? —4В-К.(5,Т,р))- ро($,Т, р)?, 
где А, В, № (57, р), К (5, Т, р) — переменные параметры, связь которых с температурой, солёностью и 
плотностью определена стандартной моделью морской воды; р(5, Т, р) =р, где Т — температура воды в 


градусах Цельсия. Скорость звука в солёной воде может быть выражена формулой: с(5.Т, р) =_/бр/бр. 


Предположим, что импульс, солёность и теплота объёма воды, ограниченного цилиндрической 
поверхностью, меняются только при изменении параметров на границе области (граничные условия Дирихле), 
в частности, изменении солёности на поверхности воды. Поведение механических параметров — плотности, 
импульса — задано условиями Неймана и внешнему воздействию недоступно. Независимо от этих допущений, 
внутри границы водной среды система уравнений (1) может быть записана в векторной форме модели 
переноса-реакции с источниковой добавкой [11, 5]: 

64 + ей =Т, 


0 дх о 


где Т= (0 Е/5+5х Ру/5+Е02Т/0х? -К(Т-Т. } — вектор, характеризующий взаимодействие потока с 


епу 


окружающей его жидкостью и планетой; а=(а ру Е) —Щ_ вектор консервативных переменных состояния; 


из <“ 
Г = (ру р-ру? у(Е+ р)) — вектор потоков, выполняющих роль обратной связи и замыкающих уравнение 


баланса. Из компонент вектора Г можно вынести общий множитель }Ё=у:4. 

На практике для решения уравнения переноса вида (6) хорошо себя зарекомендовала разностная схема, 
оператор послойного перехода которой получается линейным комбинированием аналогичных операторов 
перехода схем «кабаре» и «крест» [12—15]. Учет связи потоков массы и изменения плотности оказался более 
эффективным по регуляризованн0ой по Б. Н Четверушкину — квазигидродинамичекой — системе, 
аппроксимированной по схеме ВВЦП (разности со сдвигом вперёд и центральные по пространству). 

Дискретизацию модели сплошной среды будем проводить интегро-интерполяционным методом на 


равномерной сетке 5$ =5,х5, 5, ={х, =Ш,1=0-+и,п.й=[}, 5, ={И =, ]=0+т,т-т=Т}, где В — шаг 
сетки по вертикали, { — индекс узла (контрольного/конечного объёма в терминах метода Годунова) при 
нумерации по пространству, т — шаг сетки по времени, ] — номер временного слоя. Опишем конечно- 


разностную аппроксимацию модели (1). Разностные схемы, используемые для аппроксимации уравнений 
баланса первого порядка, дают сравнительно приемлемые результаты лишь при очень малом шаге сетки [3], что 
приводит к активному потреблению ресурсов вычислительных устройств и разработчиков, создающих 
алгоритмы и программы. 

Проблема постановки граничных условий и их согласованного задания при знакопеременной скорости 
решалась методом заполненности контрольных ячеек [7]. Учет максимальной заполненности и выражение её 
функциональной зависимостью от номера узла позволяет повысить точность аппроксимации граничных условий. 

Работа по теоретическому и экспериментальному подбору методов и схем аппроксимации привела к выбору 
методов расщепления [13] и регуляризации. Уравнения неразрывности  регуляризируется по 
Б.Н. Четверушкину. Данное уравнение подвергается разностной аппроксимации по схеме ВВЦИП. Уравнения 
баланса (переноса) импульса, солёности и полной энергии заменяются явными уравнениями, полученными 
линейным комбинированием различных аппроксимаций оператора переноса («чехарда» и «кабаре») [12, 15]. В 
соответствии с методом расщепления, перенос импульса (и скорости), солёности, полной энергии 


Информатика, вычислительная техника и управление 


217 


Бир://уезииК-аопзв.га 


218 


'А4уапсей Епртеетия Везеагсй (Коз1оу-оп-Ооп). 2023;23(2):212—224. е155М№ 2687—1653 


аппроксимирован на дробном шаге. Изменение скорости движения воды, определяющей интенсивность 
переноса массы вещества, вводится на втором дробном шаге решением волнового уравнения. Обозначим 


р=р,, р=р"1°, р=р"Н,у=у, У=у"°, у=ум, 0<о<1. В расщеплённом виде [16], после введения в 


уравнение неразрывности регуляризатора и полудискретизации на сетке 5. аппроксимации частной 


производной по времени схемой «разности вперёд по времени», первые два уравнения системы (1) могут быть 
записаны в виде: 


бу-ру 9[ д» 
РА, (рь), =-ра$ (== |, 7 
о (ру), =-р85 (х) мрт (7) 
1 р, —р, а _ Р"н —2Р» + Р"" ь _ 
=“ +(69), |= Р5 Е ‚р=Л(рТ),с=с(5Т,р), (8) 
РР р. (9) 
т 


Уравнения (7) и (9) представляют собой дискретный аналог модели переноса и изменения импульса — 
второе уравнение системы (1). Уравнение (8) прогнозирует изменение поля давления при учёте неразрывности 
потока и информации о межузловых (межъячеечных) и граничных потоках. При допущении постоянности по 
времени на каждом шаге скорости звука уравнение (8) может аппроксимироваться неявной разностной схемой 
по времени и решаться прогонкой. Именно это направление в численном моделировании и было выбрано для 
исследования. 

Уравнение (7) и четвёртое уравнение (1) аппроксимируем на эквидистантной сетке явной разностной 
схемой [13]. Слагаемые (8), определяющие волновые свойства при распространении импульса (правая часть 
уравнения), аппроксимируем центральными разностями по пространству: 


1 —1 1 1 1 1 
Код рен — 2ре + р’ х Р-Р" у Р-Р (10) 
с? д ы 2 2. 2 ь 
где К, —_ заполненность ячейки, находящейся справа от узла с индексом а К, —_ заполненность ячейки, 
находящейся слева от узла с индексом 7 - К —_ степень заполненности контрольной области Х: 1/2 < х < Х; 1/2 , 


находящейся в окрестности узла с индексом 1 (К, ‚ К, ‚ К, — переменные величины, зависящие от индекса 1) [7]. 


Оставшиеся два слагаемых (8) определяют связь скорости изменения плотности с её потоком. В (8) входит 
отношение скорости изменения плотности к отрезку времени т, выражаемой разностью со сдвигом вперёд: 


+1 _ 
р-р 


0,1 2 7 


[* 


поток массы переносится в правую часть (8) и аппроксимируется центральными разностями: 


19(58) | бор -Ф),, (ри) - Фа), 


0,1 х ах _ Вы Ъ.х К. т (11) 
где К, ,К — степень заполненности областей, находящихся в окрестности ячеек с номером 1[7]; К, 
характеризует заполненность области [х1,х..|; К — [хьхи|, № — [хх]. 
Уравнения баланса механического импульса и давления (9) аппроксимируем ВВЦП: 
к Роу Р-Р", РР 2) 


о в а 
: т ‚ 26 26 
Представим задачу разрешения разностной аппроксимации уравнения (8) в форме задачи матричной 
прогонки с переменным во времени вектором в правой части Ах=ЁР: 
Ах, +Сх+Вхи=Е. 


Аппроксимация уравнения (8), определённая на трёхточечным разностном шаблоне, в форме линейной 


системы уравнений Ар”+! = {(р", р"-", (27) .) ‚ разрешаемой относительно давления (р), имеет вид: 
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Подобная ленточная матрица получается в ходе аналогичной аппроксимации двух других слагаемых (8). 

В модели могут быть учтены динамические изменения потока на границах и других параметров, таких как 
температура, солёность, кислородосодержание. Поэтому моделирование движения воды нужно выполнять в 
цикле по времени послойно с удержанием оперативной информации о параметрах как минимум на двух 
соседних по времени слоях решения сеточного уравнения. Алгоритм вычисления гидродинамических 
параметров на двухиндексной сетке по пространству и времени включает в себя: 

— построение прогноза изменения импульса по первому уравнению (7); 

— приближённое вычисление функции пространственного распределения давления как функции плотности и 
температуры с предыдущего временного слоя; 

— оценку изменения плотности по второму уравнению (8); 

—вычисление градиента давления по новым значениям плотности и температуры, корректировка 
распределения импульса по третьему уравнению (9); 

— нахождение нового распределения полной энергии и температуры по аппроксимированному третьему 
уравнению системы (1). 

При алгоритмизации метода решения системы расщеплённых уравнений введены: 

— двоично-числовые маски, предопределяющие переключение шаблона разностной схемы с изменением 
знака скорости; 

— переменные сдвиги индексов соседних узлов, позволяющие записать решение уравнения (12) для 
граничных и внутренних узлов единой системой вычислительных операций. 

Такие переменные сдвиги индексов используются при вычислениях аппроксимаций градиентов как в 
объёме модели сплошной среды, так и на границах со вторым порядком точности по дискретным аналогам 
уравнений (7)—(9), аппроксимированных по ВВЦП и по линейной комбинации схем «чехарда» и «кабаре» 
(САВАВЕТ) [12—14]: 


47 -а к —а"'! + Ч + а 


п+1 — (п 4 п — п-1 
Е - [2 - му т - Е 0,у 20, 
а 4" 4 9" а 4" > 9" 97 о И В Ч! (14) 
1 7 3 +1 р ул 141 7 4 1 1 у! 141 1-Е — 0, ул < 0, 
т 3 2 й 3* ЗА 


где характеристика 4 для переноса соли (а = 5 ) и импульса д=ри. 


Параметр ш представляет собой «переключатель» потока при изменении знака скорости в разностных 
аппроксимациях градиентов давления и потока массы: 


0, у< 0, 


ттт = 50 


При использовании такого переключателя аппроксимация градиента при решении задачи переноса для 
импульса и солёности может быть записана формулой: 


04 х 2 И: — Ч и 4—4; +1 Ч Ч , 


Ох 3 й й 3 2й 


х=х, 


(15) 
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Применение метода заполненности ячеек позволяет корректно аппроксимировать граничные условия при 
размещении параметров состояния воды в памяти компьютера массивом числовых значений, а переменные 
сдвиги индексов позволяют уменьшить объём символьной записи подпрограммы выполнения послойного 
итерационного изменения переменных состояния, отнесённых к узлам разностной схемы [17]. 

Матричное уравнение с ленточной трёхдиагональной матрицей А (13) решается прогонкой [18]. Для 
начальной верификации программной реализации использовалась упрощённая модель состояния солёной воды 
в виде уравнения П. С. Линейкина. 

Результаты исследования. Моделирование выполнялось на основе программного комплекса, написанного 
на языке МаНаб. Отладка проводилась с использованием интерпретатора СМО Осауе. Эксплуатация также 
предполагает наличие библиотек этой системы. Программный комплекс состоит из 16-ти функциональных 
модулей. Выбор этого интерпретатора и соответствующего языка обусловлен возможностью записать и 
проверить программу, оперирующую с массивом переменных состояния, являющегося проекцией искомых 
функций р(х,), у(х,1), &(х,0) на пространственную сетку. На этом этапе исследования авторы 


абстрагировались от специфики выполнения поэлементных операций над массивами действительных чисел. 

Программная система состоит из взаимосвязанных подпрограмм: 

— выполнение одного шага по уравнению переноса по формулам (10); 

— расчет плотности морской воды (ипезсо_лгз); 

— расчет давления морской воды, находящейся в поле силы тяжести, при заданной температуре, солёности и 
плотности по эмпирическим уравнениям состояния (оТ$2Р); 

— вычисление скорости звука в зависимости от температуры, солености и давления по стандарту и 
уточнённой формуле «ЮНЕСКО» (зрее4_оЁ_зоипа); 

— оценка величины силы вязкости жидкости при движении её объёмов под действием сил давления во всех 
точках пространственной стеки (ЕогсеОЁЕйсйоп); 

— циклическое варьирование переменных состояния и времени при вертикальном движении солёной воды 
(адиа_ргосез$); 

— начальная установка константных величин, характеризующих жидкость, начальные и граничные условия 
её глобального гидрофизического равновесия (бат, зе _рагатеегз); 

— формирование ленточной матрицы, аппроксимирующей уравнение, содержащее давление (Ёапс5); 

— решение матричного уравнения методом прогонки (гип_5\еер_зНиШе); 

— преобразование температура <> полная энергия (ТЕготЕ, ЕЕготТ); 

— расчет баланса полной энергии (То{а!Ро\ег); 

— решение задачи диффузии-конвекции-реакции, в том числе при аппроксимации комбинации разностных 
схем «чехарда» и «кабаре» (АБК _зоуег); 

— оценка величин градиентов и производных по времени по центральным и направленным разностям с 
учётом изменения знака скорости и шаблона разностной схемы (41123). 

Построенная модель использована в тестовом запуске для оценки изменения плотности при повышении 
солёности поверхности водной среды на 1% на 30-Й секунде от начала моделирования и максимальной 
глубине 1км. Переходные процессы по плотности во множестве отстоящих друг от друга сечений 
(функция р(Е, х), 5, ), вызванных мгновенным изменением концентрации соли, конечны во времени (рис. 1), а 


фронт импульса плотности сильно ослабляется при движении в пространстве. 
р, кг/м 


120 130 


1010 


Рис. 1. График зависимости плотности от времени при резком изменении солёности на поверхности воды 
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Обсуждение и заключение. Представленная пространственно-распределённая модель пока не позволяет 
прогнозировать устойчивое перераспределение плотности и изменение градиента солёности при численных 
расчетах, потому что в ней не учтена плавучесть менее солёной воды и еб изменение при осолонении верхних 
слоёв. Математическое и программное обеспечение, значительно упрощающее моделирование процессов, 
приводящих к наблюдаемым эффектам галоклина, хемоклина и пикоклина, проверено и отлажено на множестве 
тестовых задач (транспорт импульса, соли, распространение Волны давления). Для высокоточного 
моделирования наблюдаемых физико-химических явлений в морях понадобится решать дополнительные 
задачи идентификации параметров модели, взяв в качестве исходной информации данные наблюдений и 
дистанционного зондирования Земли. 
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